#include<sys/stat.h>

namespace Dune {
  namespace Compose {
    template<typename EQ, typename NUMFLUX, typename INIT, typename BOUN, typename GV, typename FEMDG>
    void driver(EQ& eq, NUMFLUX& numflux, INIT& init, BOUN& boun, const GV& gv, 
      const FEMDG& femdg, Dune::ParameterTree& ptree){
      using RF = typename NUMFLUX::RF;
      static constexpr int dim = NUMFLUX::dim;
      static constexpr int nvar = NUMFLUX::nvar;
      // DEBUG
      std::cout<< " --- driver is the core part --- " <<std::endl;
      std::cout<< "dim: " << dim <<std::endl;
      std::cout<< "nvar: " << nvar <<std::endl;

      /**
       * user functions
       */
      //-- set initial condition --
      typedef InitialWrapper<dim,nvar>  InitWrapper;
      std::shared_ptr<InitWrapper> initwrapper =
        std::shared_ptr<InitWrapper>(new InitWrapper(init));
      auto uinilambda = [&](const auto& i, const auto& x) {
        return initwrapper->initial(i,x);
      };
      auto uini = Dune::PDELab::
        makeGridFunctionFromCallable(gv,uinilambda);
      std::cout<< "set initial condition" <<std::endl;

      //-- set boundary condition --
      auto blambda = [](const auto& x){return true;};
      auto b = Dune::PDELab::
        makeBoundaryConditionFromCallable(gv,blambda);
      std::cout<< "set boundary condition" <<std::endl;
      // test
      auto glambda = [](const auto& x){
        RF s=0.0; 
        for (std::size_t i=0; i<x.size(); i++) s+=x[i]*x[i];
        return s;
      };
      auto g = Dune::PDELab::
        makeGridFunctionFromCallable(gv,glambda);
      auto ascalarlambda = [] (const auto& x) {
        RF s=100.0;
        return s;
      };
      auto ascalar = Dune::PDELab::
        makeGridFunctionFromCallable(gv,ascalarlambda);

      /**
       * make grid function space for the system
       */
      //-- make grid function space --
      using VBE0 = Dune::PDELab::ISTL::VectorBackend<>;
      using CON = Dune::PDELab::OverlappingConformingDirichletConstraints;
      // using CON = Dune::PDELab::ConformingDirichletConstraints;
      // using CON = Dune::PDELab::NoConstraints;
      using GFS0 = Dune::PDELab::GridFunctionSpace<GV,FEMDG,CON,VBE0>;
      GFS0 gfs0(gv,femdg); 
      //-- make gfs for the system --
      using VBE = Dune::PDELab::ISTL::VectorBackend<
        Dune::PDELab::ISTL::Blocking::fixed>;
      using OrderingTag = Dune::PDELab::EntityBlockedOrderingTag;
      using GFS = Dune::PDELab::PowerGridFunctionSpace<
        GFS0,nvar,VBE,OrderingTag>;
      GFS gfs(gfs0);
      std::cout<< "make gfs for the system" <<std::endl;
      //-- assemble constraints --
      typedef typename GFS::template ConstraintsContainer<RF>::Type CC;
      CC cc;
      gfs.update();
      std::cout << "constrained dofs=" << cc.size() << " of "
        << gfs.globalSize() << std::endl; 

      /**
       *  make instationary grid operator
       */
      // local operator
      using LOP = Dune::Compose::DGHyperbolicSpatialOperator<NUMFLUX,FEMDG>;
      LOP lop(numflux);
      using TLOP = Dune::Compose::DGHyperbolicTemporalOperator<NUMFLUX,FEMDG>;
      TLOP tlop(numflux);
      using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
      int degree = ptree.get("fem.degree",(int)1);
      // MBE mbe(std::pow(1+nvar*degree,dim));
      MBE mbe(2*dim+1);
      // global operator
      using GO0 = Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,CC,CC>;
      GO0 go0(gfs,cc,gfs,cc,lop,mbe);
      using GO1 = Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,RF,RF,RF,CC,CC>;
      GO1 go1(gfs,cc,gfs,cc,tlop,mbe);
      using IGO = Dune::PDELab::OneStepGridOperator<GO0,GO1,false>;
      IGO igo(go0,go1);
      std::cout<< "local and global operator" <<std::endl;

      /**
       * select and prepare time-stepping scheme
       */
      // time-stepping order
      int torder = ptree.get("fem.torder",(int)1);
      Dune::PDELab::ExplicitEulerParameter<RF> method1;
      Dune::PDELab::HeunParameter<RF> method2;
      Dune::PDELab::Shu3Parameter<RF> method3;
      Dune::PDELab::RK4Parameter<RF> method4;
      Dune::PDELab::TimeSteppingParameterInterface<RF> *method;
      if (torder==1) {
        method=&method1;
        std::cout << "setting explicit Euler" << std::endl;}
      if (torder==2) {
        method=&method2;
        std::cout << "setting Heun" << std::endl;}
      if (torder==3) {
        method=&method3;
        std::cout << "setting Shu 3" << std::endl;}
      if (torder==4) {
        method=&method4;
        std::cout << "setting RK4" << std::endl;}
      if (torder<1||torder>4)
        std::cout << "torder should be in [1,4]" << std::endl;
      // set time-stepping schme
      igo.setMethod(*method);

      /**
       * coefficient vector
       */
      //-- coefficient vector, discrete grid function  CON;-
      using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
      Z z(gfs);
      using ZVDGF = Dune::PDELab::VectorDiscreteGridFunction<GFS,Z>;
      ZVDGF zvdgf(gfs,z);
      // typedef typename IGO::Traits::Domain V;
      // using V = Dune::PDELab::Backend::Vector<GFS, RF>;
      // V z(gfs,0.0);
      Dune::PDELab::interpolate(uini,gfs,z);

      Dune::PDELab::constraints(b,gfs,cc); // assemble constraints
      // Dune::PDELab::set_nonconstrained_dofs(cc,0.0,z);
      // Dune::PDELab::set_constrained_dofs(cc,0.0,z); // set zero Dirichlet boundary conditions


      /**
       * make vector consistent NEW IN PARALLEL  
       */
      Dune::PDELab::ISTL::ParallelHelper<GFS> helper(gfs);
      helper.maskForeignDOFs(z);
      // Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,z);
      Dune::PDELab::AddDataHandle<GFS,Z> adddh(gfs,z);
      if (gfs.gridView().comm().size()>1)
        gfs.gridView().communicate(adddh,
          Dune::InteriorBorder_All_Interface,
          Dune::ForwardCommunication);

      /**
       * prepare VTK writer and write first output 
       */ 
      // create directory
      int subsampling=ptree.get("output.subsampling",(int)1);
      using VTKWRITER = Dune::SubsamplingVTKWriter<GV>;
      VTKWRITER vtkwriter(gv,Dune::refinementIntervals(subsampling));
      std::string filename=ptree.get("output.filename","output");
      struct stat st;
      if( stat( filename.c_str(), &st ) != 0 ) {
        int stat = 0;
        stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
        if( stat != 0 && stat != -1)
          std::cout << "Error: Cannot create directory "
            << filename << std::endl;
      }
      // add auxiliary variables to vtkwriter
      using Vt = Dune::PDELab::Backend::Vector<GFS0, RF>;
      typedef Dune::PDELab::DiscreteGridFunction<GFS0,Vt> DGF;
      Vt x(gfs0);
      DGF dgf(gfs0,x);
      Vt factor(gfs0);
      DGF dgffactor(gfs0,factor);
      // typedef Dune::PDELab::MinusGridFunctionAdapter< decltype(g),decltype(ascalar)> MinusGFAdapter;
      // MinusGFAdapter aminusGFadapter(g,ascalar);
      typedef Dune::PDELab::MinusGridFunctionAdapter< DGF,DGF> MinusGFAdapter;
      MinusGFAdapter aminusgf(dgf,dgffactor);
      typedef FastEigenvalueDGF<GFS,Z> FastEvDGF;
      FastEvDGF fastev(gfs,z,eq);
      // vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<DGF> >(dgf,"test"));
      // vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MinusGFAdapter> >(aminusgf,"test"));
      vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<FastEvDGF> >(fastev,"test"));
      // create vtk sequence writer 
      using VTKSEQUENCEWRITER = Dune::VTKSequenceWriter<GV>;
      VTKSEQUENCEWRITER vtkSequenceWriter(std::make_shared<
        VTKWRITER>(vtkwriter),filename,filename,"");
      // add data for all components of the space to the writer
      std::cout << "output initial condition" << std::endl;
      Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter,gfs,z,
        Dune::PDELab::vtk::DefaultFunctionNameGenerator("u"));
      // V testvar(gfs,0.0);
      // typedef Dune::PDELab::DiscreteGridFunction<GFS,V> TestDGF;
      // TestDGF testvardgf(gfs,z);
      // typedef Dune::PDELab::VTKGridFunctionAdapter<TestDGF> VTKTestDGF;
      // vtkSequenceWriter.addVertexData(std::shared_ptr<VTKTestDGF>(new
      //                                   VTKTestDGF(testvardgf,"testvar")));
      // output initial conditon
      vtkSequenceWriter.write(0.0,Dune::VTK::appendedraw);
      Dune::PDELab::interpolate(g,gfs0,x);
      Dune::PDELab::interpolate(ascalar,gfs0,factor);

      /**
       * Linear problem solver backend
       */
      using LS =  Dune::PDELab::ISTLBackend_CG_AMG_SSOR<IGO>;
      int verbose=0;
      if (gfs.gridView().comm().rank()==0) verbose=1;
      LS ls(gfs,100,verbose);
      // time-stepping
      // typedef typename IGO::Traits::Domain V;
      // Dune::PDELab::ExplicitOneStepMethod<RF,IGO,LS,V,V> osm(method1,igo,ls);
      Dune::PDELab::ExplicitOneStepMethod<RF,IGO,LS,Z,Z> osm(method1,igo,ls);
      osm.setVerbosityLevel(5*verbose);

      /**
       * time loop
       */
      RF time = 0.0;
      RF T = ptree.get("initial.T",(RF)1.0);
      RF dt = ptree.get("fem.dt",(RF)0.1);
      const int every = ptree.get<int>("output.every");
      int counter = 0;
      // V znew(gfs);
      Z znew(gfs);
      double dtfix = 0.0;
      dtfix = ptree.get<double>("fem.cfl") / dim * 
        (ptree.get<double>("grid.structured.LX") / 
        ptree.get<int>("grid.structured.NX"));
      std::cout << "begin time loop"  << std::endl;
      while (time < T) {
        std::cout << "a new loop begin, loop counter: " << counter << std::endl;
        RF maxfastev=0;
        maxfasteigenvalue(gv,fastev,maxfastev); 
        std::cout << "maxfastev: " << maxfastev << std::endl;
        dt = dtfix / maxfastev;
        // getdt(maxfastev,dt,ptree);
        numflux.set_dt(dt);
        std::cout << "timestep dt: " << dt << ", current time: " << time << std::endl;
        // do time step
        osm.apply(time,dt,z,znew);
        // accept time step
        z = znew;
        time += dt;
        counter++;
        // output to VTK file every n-th timestep
        std::cout << "output snapshot " << counter << std::endl;
        if(counter % every == 0)
          vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
        std::cout << "this loop end" << std::endl;
      } // time loop
    } // driver 
  } //Compose
} //Dune
